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ABSTRACT 


The equations of horizontal motion of the neutral atmosphere be- 
tween 120 and 500 km are integrated with the inclusion of all the non- 
linear terms of the convective derivative and the viscous forces due to 
vertical and horizontal velocity gradients. Empirical models of the 
distribution of neutral and charged particles are assumed to be known. 
The model of velocities developed is a steady state model. In part 1 
the mathematical method used in the integration of the Navier-Stokes 
equations is described and the various forces are analysed. 
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FULL NON-LINEAR TREATMENT OF THE 
GLOBAL THERMOSPHERIC WIND SYSTEM. 


Introduction 

The wind system of the thermosphere has been a subject of a considerable 
number of observational and theoretical investigations during the last decade. 

A knowledge of this wind system is required for an understanding of the struc- 
ture of the neutral atmosphere and its energy balance and the ionosphere. It 
has been suggested that the solution of the phase problem (the phases of the 
density and temperature in the thermosphere) of the neutral atmosphere is in- 
timately connected with the wind system. Many ionospheric effects concerning 
the latitudinal distribution of charged particles, the daily variation of electron 
density and the maintenance of the nighttime ionosphere cannot be fully explained 
without taking into account the effect of the atmospheric wind system. 

Several computations of the thermospheric wind pattern have been made. 

All of them solve the horizontal equations of motion and assume hydrostatic 
equilibrium. There are two general approaches: (1) a perturbation treatment 
of the full set of hydrodynamic equations (Lindzen (1970); Volland and Mayr 
(1970 and 1972); (2) a solution of the equation of motion using a given model of 
atmospheric structure (Geisler (1967); Bailey et al. (1969); Challinor (1970); 

Cho and Yeh (1970 Y; Ruester and Dudeney, (1972)). The various computations 
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differ in the assumptions regarding the treatment of the ion and viscous drag. 

All of them have in common that the non-linear terms of the convective deriva- 
tive of the flow velocity were not fully included. Especially the non-linear terms 
due to meridional velocity gradients have not been included in the computation 
by any of the investigators. 

We have chosen the second approach, i.e„, we have used in our calculations 
the atmospheric model of Jacchia (1965) and the Penn State Ionospheric Model 
of Nisbet (1970) and solved the full Navier-Stokes equations including all non- 
linear terms and the full expression for the viscous forces. None of the earlier 
papers have included in the computations of the viscous terms the horizontal 
velocity gradients. 

1. MATHEMATICAL METHOD 

The Steady State 

The neutral air motions treated here are steady state motions. The steady 
state refers to the sun's system. All variables like density, temperature, veloci- 
ties and the various forces are assumed to be independent of universal time in a 
coordinate system fixed with respect to the sun, if we limit ourselves to periods 
of a few days. Variations with universal time due to the changing declination of 
the sun, solar activity or the semi-annual density variation are not excluded, 
but each change of this type is treated as a separate model calculation. There- 
fore for each day of the year and a given solar activity our set of equations has 


2 



5 


three independent variables. For these we have chosen the height z, i.e. the 
distance from the surface of the earth, the geographical co- latitude 6 and the 
azimuth 4> referred to a coordinate system not rotating with the earth. 4> is 
identical with the local time r. These simplifying assumptions exclude the 
treatment of physical processes that have a true universal (not local) time 
dependence with time scales of a day or less. Such processes would include all 
effects that are due to the inclination between the earth’s axis and its magnetic 
dipole. In the reference frame that rotates with the earth the independence from 
universal time means that all variables depend only upon the combination k-o>t 
of geographical longitude k and universal time t, i.e. all phenomena are peri- 
odic with the period 2 tt/co of one solar day. Non-recurring changes, sporadic 
variations and variations that depend explicitly on longitude k and not only on 
t = k-co t are excluded from our treatment. 

To obtain the appropriate form of the hydrodynamic equations for an observer 
on the earth, we must first transform the equations to a frame rotating with the 
earth and then take the limit to the steady state. Let 4> be the longitude in the 
fixed frame and X the longitude in the rotating frame, and let t' be the time in the 
rotating frame and t in the fixed frame. Then the transformation equations are 

k = <f> - cot (1) 

t' = t 

The partial derivative with respect to 4> and t are according to the chain rule 
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_ 3 _ _ _ 3 _ 

3 <£ 3 \ 

_3 _ 3 3 

3 1 3 1 ' 


( 2 ) 


In the steady state we can let the partial derivative with respect to the time t 
in the fixed frame go to zero and thus obtain the following relationship between 
the time variations and the longitude variations for an observer on the earth 

(3) 


3 _ 3 _ 

3 1 ' 


Care must be exercised when this relationship is used; in particular one must 
not try to describe transient phenomena using equations based upon these rela- 
tionships. Boundary conditions imposed upon the equations must be consistent 
with the assumption of the steady state. 

It is to be noted that the solutions of the tidal equations (Chapman and 
Lindzen, (1970)) that would correspond to the steady state are those for which 
the parameters f and s (in their notation) of tidal theory have the ratio unity. 


Equations of Motion 

The Navier-Stokes equations in vector form in the rotating frame are 
D v 

+ 2 a) x V - rj div grad V - r//3 grad div V + f . Qn = - (grad P )/ p + g (4) 

where V is the velocity vector, S the rotational velocity of the earth, f ion the 
ion drag force, P the pressure, p the density and 77 the kinematic viscosity. 

We are concerned with the global wind pattern with periods of the order of 
one day and its harmonics. The appropriate coordinate system for this problem 
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is spherical. As the vertical velocities are smaller by more than one order of 

magnitude than the horizontal velocities, we may solve the horizontal equations 
neglecting the vertical velocity terms that appear in the convective derivative, 

the Coriolis force and the viscous forces. 

Let be the meridional velocity, measured positive from North to South, 

and V ( ^> the zonal velocity, measured positive from West to East, then the 

convective derivative of the velocities becomes 


DV< 5 > 3V< fl > 1 

= co + — 

Dt 3 cjo r 



30 


V(^) 3V (i9) 
sin 0 3 cp 


( V C^))2 ctn g\ 


DV^ 3V^> 1 

= co + — 

Dt • dcfi r 


J(.9) 


3V<*) _VW 3V<*> y(d) yofi) ctn 0 

1 • r\ / 1 


30 


sin 0 3 <b 


(5) 


with 0 the co-latitude and 4) the longitude in the sun's system (or local time in 
the earth's system) and r the distance from the center of the earth. 

We have made use in these equations of the steady-state assumptions, i. e., 
we have assumed all variables to depend only upon local time and not explicitly 
upon geographic longitude and universal time. 

Substituting the expressions for the convective derivatives in the Navier- 
Stokes equations and dividing byo>, it is seen that the non-linear terms are 
multiplied by a factor l/(«r). The numerical value of cor is 462 m/sec. The 
main inertial term 3V/3^is of the order of the velocity itself. Without detailed 
calculations this suggests that for velocities that are lower by an order of 
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magnitude than ^r it is to be expected that the non-linear terms have no marked 
influence on the results,, On the other hand for velocities that are comparable 
or larger thanwr it is essential to include the non-linear terms. This estimate 
is based on the assumption that the derivatives of the velocities are of the order 
of the velocities themselves (in the spherical coordinate system) which is ob- 
viously true for the zonal derivatives, but may not be true, and indeed is not in 
the equatorial region, for the meridional gradients. 

The equations we now have to solve are a pair of non-linear coupled second 
order partial differential equations with the three independent variables r, Q 
and 4> or r, 9 and r. According to the steady state picture the solutions must be 
periodic in local time (the variable r), finite and continuous everywhere. The 
boundary conditions in altitude specify a given distribution of velocities at the 
lower boundary and a zero vertical gradient of horizontal velocities at the upper 
boundary. Not every density model and ion distribution would result in a steady 
state solution. This is seen by the simple considerations that pressure gradients 
that are too large cannot be maintained due to the rapid flow which would tend 
to equalize these pressure gradients. We shall assume that a solution exists 
based upon empirical models of the atmosphere and ion distribution. From the 
assumption of the existence of the solution certain limitations of the velocity 
distribution may be obtained. At equinox conditions, where from symmetry 
considerations the meridional component of velocity is zero at the equator, the 
equation for the zonal component of velocity becomes 
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Bv 

3 T 


1 Bp B 2 V 

D (r) V + 7] (r) 


i r p B r 


( 6 ) 


where D is the ion drag coefficient for unit mass. 

A theorem by Kallina (1970) states that a sufficient condition for the existence 
of a periodic solution of equation (6) is that the coefficient (1 + v/wr) of 
Bv/B r is positive for all t and r. Thus, in this simple case, we may expect a 
periodic solution with the model pressure gradients and ion drag coefficients 
as long as the zonal velocity in the westward direction remains less than the 
earth’s rotational velocity. However, if the zonal velocity exceeds the rotational 
velocity of the earth in the westward direction then one would expect to obtain 
exponential growing solutions and no periodic or steady state solution. A steady 
state solution for such a model of the atmosphere and ionosphere would not 
exist. If at any stage in the computational process one obtains zonal velocities 
of this behaviour numerical instabilities will develop. For example, if the values 
of the driving force divided by &> would be used as initial values for V in an itera- 
tive procedure for the solution of the steady state horizontal flow equations, 
then with the model pressure gradients of the Jacchia model the numerical 
solution would diverge. Kallina’s theorem needs to be generalized in order to 
guide us concerning the existence of a solution for the full set of equations (4) 
for the horizontal flow under general conditions. Such a generalization is at 
present not available. 
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In a purely linear treatment this property of the velocity field would not be 
apparent. Any model of atmospheric density and ion distribution results in a 
steady state velocity field when the non-linear terms of the convective derivative 
are neglected. Only the full non-linear treatment shows that the steady state 
does not exist if certain conditions are not met, i. e. the ratio of the driving 
forces to the drag forces is too large. It may be argued that Kallina’s theorem 
does not unequivocally show this property of the velocity field, as it treats only 
a very simplified system of equations and furthermore gives only sufficient and 
not necessary conditions for the existence of a steady state solution. On the 
other hand we have made a large number of numerical model calculation that 
indicated that a generalization of Kallina’s theorem must be true. These model 
calculations converged always when a parameter R, which indicated the ratio 
between driving and drag forces, was smaller than a certain R 0 . They diverged 
for all models with R greater than R 0 , thus indicating that steady state solutions 
for such models do not exist. The value of R 0 was well defined, a very small 
excess of R over R 0 caused the models to diverge immediately. 

Method of Solution 

The condition of periodicity will be satisfied if we expand each term of 
equation (4) into Fourier modes with respect to local time. The velocity is 
decomposed into average (or zero order), diurnal and semi-diurnal harmonic 
components; each component is taken as a function of latitude and altitude. 
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Also the pressure gradients divided by density, the kinematic viscosity and the 
ion drag coefficient are expressed in a similar fashion. Accordingly the veloc- 
ities will be expressed by 


l 

V (6,) (0, 9, r) = (a£ 5) (9, r) cos 2vk4> + ( 6 , r) sin 27xk 4>) 


k a o 


(7) 


•((£, Q t r ) = y 1 (9, r) cos 277k<£ + (9, r) sin iTTkcfr) 


k = 0 


All products of Fourier terms which give rise to modes higher than the semi- 
diurnal mode are dropped in the calculations. Thereby we have reduced the 
problem to the solution of ten coupled non-linear partial differential equations 
of second order for the ten unknown functions a k and b k (with b 0 = 0) with altitude 
and latitude as the independent variables. A finite difference scheme is used to 
obtain the solutions. The derivatives with respect to latitude are replaced by 
second order differences, except at the north and south poles, where forward 
and backward differences are used respectively. The first derivative with 
respect to altitude is not present in our equations. The second derivatives with 
respect to altitude are replaced by normal second order differences. The bound- 
ary conditions in altitude are expressed to second order accuracy (Varga, p. 191). 
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Treatment of the Equations at the Poles 

Previous investigations of the global wind field have either excluded the 
polar regions of have used a simplified computational scheme (sometimes using 
Cartesian coordinates). This can be done without great loss of accuracy if the 
latitudinal coupling of the wind field is neglected. In our treatment such simpli- 
fications are not possible. The inclusion of the poles in the computational 
scheme requires some care and we therefore present some details of our 
approach. • 

The boundary condition that the solution be single valued at the poles re- 
quires that all modes of the horizontal velocity but the diurnal mode vanish at 
the poles as can be seen by simple kinematical considerations. Furthermore, 
the meridional and the zonal components of velocity are ninety degrees out of 
phase. Thus we must have 


a ( ^= b ( ^- sig 


(8) 


(«) 


'j - --! Slg 

where sig = +1 at the north pole and -1 at the south pole. Also we have 

b <S) = a<*>= b<*>= n 

k k K k 

for k not equal to one. These kinematical boundary conditions at the poles 
reduce the number of unknown velocity components at each pole from ten to two. 

On the other hand the system of equations (4) yields four separate equations 
(two for each of the diurnal modes for both the meridional and the zonal com- 
ponents), so the problem seems to be overdetermined at the poles. 
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Formally at the poles the Navier-Stokes equations appear to have a singu- 
larity due to the factor (l/sin0) that appears in some of the iion-linear terms. 
Let L (0) and I ^ be those terms then 


= 


V^) 
sin 6 


BV< g > 
B c f> 


- cos 9 


L<*> = 


sin 9 


BV^> 
B 4> 


+ cos 9 


( 10 ) 


But from boundary conditions (8) and (9) it is seen that the bracketed parts of 
L (i9) and L (<j!,) vanish at the poles and L (S) and become therefore inde- 
terminate. To evaluate this indeterminate form of the non-linear terms we 
have to apply L'Hospital’s rule at the poles. The expressions (10) become 


= 


V^) 
cos 9 


B 2 V< 0 > „BV^> 

- cos 9 — — — — 

B 9 B cp B 9 


L<*> = 


y(^) 
cos 9 


B 2 

d9d(p 


+ cos 9 


BV (5) 

B0 


( 11 ) 
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We have to evaluate L ( ”' ) and L — for the diurnal modes only, as we have 
already shown that the other modes vanish at the poles. It is easily seen that 
the product d 2 V/dO dcfiis zero also for the diurnal modes. Therefore at 
the poles the expressions (11) are reduced to 


L<«> - - v« 5V<1> 
dO 

( 12 ) 

l(^) = 

dO 


Substituting (12) into the full expression of the convective derivative as given 
by equations (5) we obtain 


DV/ 0 ) 

Dt 


dcfi 




do do J 


£!£= « ?L!!! ♦ I /U 2^ 1 1 v <*> ^ 

Dt d<p r ^ dO- dO y i 


(13) 


The subscript 1 means that the diurnal modes only are to be taken. Imposing 
the additional requirement that the derivatives with respect to of the semi-diurnal 
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components vanish and again making use of the properties of the Fourier com- 
ponents of velocity at the poles as given by (8) and (9) it is seen that the four 
expressions for the convective derivative are reduced to two expressions. All 
other terms of the Navier-Stokes equations, i. e. the driving forces, drag forces, 
Coriolis forces etc. are also vectors in the horizontal shell given by r = constant. 
For this reason the same kinematical conditions expressed by (8) and (9) hold 
for all the terms of (4). Thus, from these considerations it is seen that at each 
pole the four Navier-Stokes equations that remain (two for each of the diurnal 
modes and velocity components) are really only two independent equations. The 
over-determination of the problem at the poles is therefore only apparent and 
not real. By using the procedure outlined above we have guaranteed that our 
solutions will satisfy the Navier-Stokes equations at the poles and will be con- 
tinuous over the whole globe. 

Although the vertical velocity does not appear in our set of equations, it 
should be remarked that the behaviour of the Fourier components of the vertical 
velocity at the poles is completely different. All the modes of the vertical 
velocity except the average (or zero) mode vanish identically at the poles. 

The Iteration Procedure 

The problem has now been reduced to a system of non-linear algebraic 
equations where for each mesh point in latitude and altitude there are ten 
unknowns and ten equations. If latitude differences of five degrees and altitude 
differences of 20 km are chosen then for the altitude range of 120 km to 500 km 
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there result 7400 equations with 7400 unknowns. These are reduced by our 
boundary conditions at 120 km to 7030 unknowns. These equations are repre- 
sented symbolically by 

F(V) = R < 14 ) 

where F represents a non-linear matrix operator on the 7030 unknowns repre- 
sented by the vector V. The vector R on the right hand side of (14) is derived 
from the pressure gradients and also has 7030 components. The method adopted 
to solve the above system of equations is a double iteration scheme consisting 
of a single Newton-Ralphson procedure combined with a Gauss-Seidel iteration 
(Ortega and Rheinboldt, p. 214). To apply the Newton-Ralphson technique we 
expand about a previous iterative value V < 1 > , where for i = 0 we use some 
initial estimate (usually zero). This may be expressed as 

F(V (i) ) + (V(i + !) _ V (i) ) = R (15) 


or 

^ITz-V fil 

-V .. ' .. (VO+ 1 ) - V(O) - R _ F (Y (i) ) (16) 

j 

where 3F/BV is the Jacobian J or Frechet derivative. We shall denote yO+i) - 
V (l) by W (l + 1) c Equation (16) results in the formal solution 

y(i + l) = y(i) + J-1 (R _ F( y(i))) ( 17 ) 
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or 

W< i+1 > = J" 1 (R - F(V (i) » (18) 

Obviously the direct inversion of the matrix J or the solution by an elimination 
process is not possible with the present computational facilities. To obtain a 
solution a double iteration scheme is required. 

To accomplish this we may order the vector V and the Jacobian J in either 
of two ways. The first method (hereafter referred to as method one) is illustrated 
in figures 1 to 4. The matrix J is ordered into blocks where each block cor- 
responds to a certain altitude. The blocks themselves are ordered into sub- 
blocks according to latitude. These sub-blocks are matrices of order 10 x 10 
according to the ten Fourier modes of the velocities. With this ordering the 
Jacobian Matrix has the form of a tri-diagonal block form matrix in which the 
non-linear terms appear only in the diagonal block. A Gauss-Seidel iteration 
procedure is adopted, in which the matrix J is decomposed as follows 

J = D - L - U (19) 

where D,L, U are the diagonal, the lower tri- angular and upper tri- angular 
matrices respectively. The matrices Land U in this case contain terms arising 
only from the viscous coupling in the vertical direction. Let the index k denote 
a fixed altitude and D k , L k and U k the sub-blocks of D, L and U corresponding to 
the index k. Then the system of equations (16) can be written as 

D k (V°) W k (i + 1) =L k W k-i (i + 1) + U k W k + i (i) + (R-FCV^- 1 ))), (20) 

with k running from one to nineteen. 
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The matrices L k and ll k have diagonal block form. Equation (20) represents a 
system of 370 non-homogeneous linear algebraic equations. The right hand side 
is known explicitly. Due to the tri-diagonal form of D k the solution can be ob- 
tained by a tri-diagonal block elimination scheme (Vargas 1962, p. 196). 

The matrices D k , L k and U k are computed from algebraic expressions. 
These algebraic expressions are obtained by differentiating the system of non- 
linear algebraic difference equations with respect to the unknowns, i.e. the 
components of V. The construction of these expressions, as well as the forma- 
tion of the algebraic equations themselves, was performed with the aid of Formac 
(Sammet, 1967), a computor program which performs algebraic manipulations. 

An alternate scheme (method 2) is to order the vector V and the Jacobian 
matrix according to 37 blocks of latitude. Each of the blocks of size 190 x 190 
is ordered according to altitude. Thus the index k in equation (20) denotes a 
fixed latitude. In this scheme the matrices L k and U k contain terms arising 
from the non-linearities of the Navier Stokes equations. The matrix D k is again 
of tri-diagonal block form and the solution can be obtained by tri -diagonal block 
elimination as described above. An advantage of this method (method 2) is that 
with an initial value of zero for the velocities of the start of the iteration scheme 
the velocities obtained from the first iteration are the solutions of the Navier- 
Stokes equations when the non-linear and horizontal viscosity terms are ignored. 
This linear solution will contain the full effects of vertical viscosity. Thus it 
corresponds to the results of previous investigations (Geisler, 1967; Bailey et 
al, 1968). 
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The first method of ordering converged in all cases independent of mesh 
size, though it required a large number of iterations. The second method of 
ordering did not converge in all cases studied. Both methods yielded the same 
solution when convergence was obtained. The second method was applied to 
mesh sizes of five and ten degrees in latitude. Convergence was obtained at 
equinox conditions for a ten degree mesh without difficulty. To obtain conver- 
gence for the five degree mesh at equinox conditions a convergence factor 
(Ortega, 1970, p. 187) was added to the diagonal elements of the Jacobian. Con- 
vergence was obtained at solstice with a ten degree mesh in latitude also, by 
adding a convergence factor to the Jacobian. No convergence was obtained at 
solstice conditions with a five degree mesh in the second method, although 
various weight factors and smoothing techniques were attempted. Also the 
Jacobian method of iteration (Ortega, 1970) was attempted in this case with no 
better success. 

The second method of ordering was used to test the accuracy with respect 
to changes of the altitude mesh size and to study the effects of various ion drag 
and viscosity coefficients. 
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2. ANALYSIS OF FORCES 
The Driving Force 

The right hand side of our system of equations (4) are the driving forces 
f that arise from the horizontal pressure gradients. The driving force is given 
by 


f d = - 7 v „ p < 21 > 

We may derive f d from a given atmospheric model that specifies the density and 
temperature distribution as a function of latitude, local time t and height z., 

In the lower thermosphere the driving forces are relatively small and therefore 
the geostrophic approximation of the equations of motion given by 

2 [SxV] = --L V P (22) 

P h 

yields fairly accurate results for the wind field. The winds are controlled by the 
Coriolis force and are perpendicular to the pressure gradients as is apparent 
from equation (22). In the thermosphere this is not true as we may no longer 
neglect the other terms of the equation of motion, especially the inertial terms. 
The winds that result from an integration of the equation of horizontal motion 
are in a direction that is close to the direction of the pressure gradients, as will 
be shown by our results and also by previous work on this subject. For this 
reason the global pattern of the pressure gradients is the most important param- 
eter in the determination of the thermospheric global wind pattern. We s hall 
specify some of the characteristics of the pressure gradients used in our compu- 
tation. 
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Similar to most previous calculations (Geisler, 1967, Kohl and King, 1967) 
we have chosen the Jacchia model for the determination of the pressure gradients 
as this model yields .a density distribution that is in close agreement with satel- 
lite drag derived densities in the isothermal region. The forces are derived by 
a differentiation process from the quantities described by the atmospheric model, 
therefore the reliability of the forces so derived will be considerable less than 
that of the density and the temperature given by the model. In particular we 
would like to point out the uncertainties that arise when the Jacchia model is 
used. 

(1) In the Jacchia model the density p and temperature T are continuous 
functions of the independent variables 6, r and z. This is equally true for the 
poles. On the other hand the derivatives of density and temperature are not 
continuous at the poles and therefore the driving force remains undefined at the 
poles. For a global description of the driving forces it is therefore necessary to 
modify the Jacchia model in the polar regions. We have used such a modification 
(Blum and Harris, 1973) to overcome this difficulty. 

(2) The temperature in the Jacchia model is a model parameter and not 
necessarily identical with the true kinetic temperature. For this reason the 
pressure calculated from Jacchia’s model may deviate from the true pressure 
both in amplitude and in phase. If incoherent back scatter radar observations 
of temperatures are used in place of the Jacchia model temperatures, it may be 

V 

expected that a more realistic pressure pattern would be obtained. Such a 
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pressure pattern would introduce a phase shift in the wind pattern of approxi- 
mately half the phase difference between radar temperature and drag derived 
densities. 

(3) At present the data on the density and temperature distribution are 
inadequate to form a definite global model in the height region between 120 to 
250 km. Therefore the pressure gradients deduced from the Jacchia model for 
this height region are extremely unreliable and no great validity should be at- 
tached to the computed wind pattern below 250 km. The effect of this uncertainty 
of the driving force and with it the flow pattern, on the isothermal region is not 
very considerable. The reason for this is that the coupling between the flow 
patterns at various heights is only through the viscosity term ^ 2 y /b z 2 . This 
term is not very important in the lower thermosphere as the kinematic viscosity 
has a relatively low value in this height region. Due to its exponential increase 
with height it only becomes important at higher altitudes (Figure 12). The in- 
fluence of the flow below 250 km on the flow on the isothermal region will there- 
fore be small. The pattern in the isothermal region is essentially determined 

by the forces that exist in this height region and not by the flow at lower 
altitudes. 

(4) Relatively minor modifications of the density and temperature distribu- 
tion especially as regards their latitudinal dependence which is less well known 
that their local time dependence, may have considerable effect on the meridional 
driving force and thereby on the meridional flow pattern. Drag data do not 
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permit a determination of the latitudinal variations of the atmosphere that is 
sufficiently accurate to conclusively ascertain the meridional pressure gradients. 
A slightly different model of the latitudinal density dependence would result in a 
substantially altered meridional force. The wind pattern, especially the meri- 
dional winds, that would result from the pressure gradients based on the densities 
determined from the OGO 6 mass spectroscopic observations (Hedin et al. 1972) 
would be different in essentials from the winds derived with the driving forces 
based on Jacchia's model. 

We shall describe some of the properties of the driving forces that result 
from Jacchia's model: 

(1) At equinoxes the Jacchia model is symmetric with respect to the equa- 
tor. For this reason both amplitudes and phases of the azimuthal driving forces 
are equal for the corresponding points of the two hemispheres. The meridional 
forces for corresponding points are equal in amplitudes but have a phase differ- 
ence of 12 hours, i.e. the forces at corresponding points are both directed either 
towards the equator or the respective pole. Furthermore, the diurnally averaged 
meridional force (i.e. the Fourier component of order zero) is directed towards 
the equator for both hemispheres. It has a maximum amplitude at a latitude be- 
tween 40° and 30°. This result differs from Rishbeth's result (1972) who has 
determined a vanishing diurnal average meridional force from the Jacchia model. 
On the other hand the amplitude of the diurnal variation of the meridional force 
is always larger than the diurnally averaged component as shown by figure 5 for 
the northern hemisphere. 
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(2) At summer solstice the diurnally averaged meridional forces are 
directed towards the south pole for all latitudes. They have a maximum near a 
latitude of 30°. At winter solstice the direction is reversed and the maximum 
of the force is at -30°. In contrast to equinox conditions the amplitude of the 
diurnal variation of the meridional forces is less than the diurnal average ampli- 
tude in the latitude region between 30° to -20° for summer solstice conditions. 
Therefore the meridional forces in this latitude belt point during the whole di- 
urnal cycle to the south in summer and north in winter . Figure 6 shows this 
behaviour of the driving forces at summer solstice conditions. 

(3) In the Jacchia model the shape of the local time variation of the pressure 
and especially its extrema, are independent of latitude. From this follows im- 
mediately that all Fourier components of the azimuthal and meridional driving 
forces are 90° out of phase. 

(4) We shall investigate the dependence of the driving forces on height and 
solar activity. It will be seen that in the isothermal region a first approximation 
yields driving forces that increase almost linearly with height and are nearly 
independent of solar activity. These properties of the driving forces can easily 
be demonstrated by numerical computations based on the Jacchia model. With 
respect to the linear height dependence this is shown in figures 5 and 6. Analyti- 
cally these properties become apparent when the slight dependence of the driving 
forces on the variation of the mean molecular weight M with latitude and local 
time is neglected. In order to show the dependence of the driving forces on the 
separate variations of temperature and density we may rewrite equation (21) 
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- I V p = _ V in p - - V T + 51 V, M 
- h M h M h M a h 


P 


f . = -gHV. In P = -gHV h (1 




dz/H ) 


(23) 


(24) 


with H the mean scale height and R the universal gas constant. 

Using (24) as our starting point we note that in the Jacchia model the pres- 
sure P at the lower boundary does not depend on geographical position, local 
time and solar activity. Therefore 


= g H 


I. dz 


(25) 


where we have replaced V h by B/3x for convenience of notation and x is any 
horizontal coordinate. The forces for to a single atmospheric constituent become 


7j = eT l £(f-) d * = ‘ T £f 


dz 

IF 


(26) 


The variations of the temperature T are derived from the Jacchia model 


T( z > x) = ^(x) - (T ro (x) - T 120 ) exp (- cr (z - z Q » (27) 

where a is related to the temperature gradient at the lower boundary and is only 
slightly dependent on To,, the exospheric temperature. To, itself is dependent on 
x (local time and latitude) and on the solar activity S 

T m (S.x) = T 0 (1 + aS) r(x) (28) 
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where a = 3.24 (flux units) -1 and not depending on either x or S, T 0 = 383° an 
extrapolated value of the exospheric temperature for S = 0 and r(x) describes 
the dependence of the exospheric temperature on latitude and local time. As will 
be seen in the following, the near independence of the driving forces on solar 
activity is due to the possibility to express the exospheric temperature by eq. 
(28) where the function r is independent of S. Equations (27) and (28) make it 

possible to evaluate the integral dz/T in closed form with the result 

z o 

| dz/T= (z - z 0 )/T m (S,x) + In ^Ii^.j/(T co (S > x)a) (29) 

The force f d becomes for the isothermal region, where T is replaced by T m , 

7 d = - B'/X, ' ^ • - 3 -(ln (30) 

The second term of (30) is independent of the height z. The linear dependence 
of the forces on altitudes results immediately from (30) with the coefficient 
~g/T<x> 3T (Jj /ox . The variation of the forces with solar activity is given by 

3f d /3S = gz in T ro ) - g/cr JL T*, 1 - (In (T ro /T 120 )/TJ (31) 

The first term of (31) vanishes due to the form of T^ according to (28). The 
second term must be calculated explicitly: 
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si (‘"fej-ar 11 " 


( T \ 3T 

V 1 /T ”) (32) 

V T 1 20/ Bx 


The first term of (32) vanishes again because it is a mixed derivative of In . 
A further evaluation of the second term of (32) yields the final result 


3? ,/3S = 1 


/ 

r a 
r 1 + aS 


(33) 


In order to evaluate the relative change of f d with solar activity we use the 
simple estimate 



derived from (26) and therefore 


(34) 



3S 


1 a 

crz 1 + aS 


(35) 


For the determination of the relative change of f d given by (35)'we have to use the 
constants of the Jacchia model {a = 0.03 km" 1 ). For a solar activity S = 200 
and a height z = 300 km we obtain 


1 


111 

3S 


1_ 

1800 


(36) 


A change of 50 flux units would therefore cause only a fractional change of the 
forces of about 3% which is insignificant. While our derivation is not strictly 
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correct for the real atmosphere which has several constituents, it is a good ap- 
proximation for the isothermal region because there the driving forces are 
mainly due to the variations of atomic oxygen which is the major constituent in 
the isothermal region below 500 km. The exact numerical computation for the 
dependence of the forces on solar activity bears out the above estimates: below 
250 km the dependence of the driving force on solar activity is complicated and 
no clear trend is easily discernable, above 250 km the near independence of the 
driving force on solar activity is reaffirmed by the numerical results. 

The local time and latitudinal distribution of the driving force at a height of 
300 km is shown for both equinox and solstice conditions in figures 7 and 8. 

The Ion Drag Force 


The ion drag force f . on per unit mass is given by 


f. 


N. 


v. (V. - V') L 

in ' ion v / n 


(37) 


— > 

where V. on is the ion velocity and v in the ion- neutral collision frequency for 
momentum transfer. The collision frequency is proportional to the density 
and depends slightly on the types of ions present (Stubbe, 1968). The drag force 
is therefore determined by the ion velocity and the ion distribution. We have not 
used the individual distributions for each ion species for the determination of 
f ion .butused the approximate expression (Chapman, 1965) that gives the ion- 
neutral collision frequency as a function of N the total density. This 
expression is 


i^ in /N = 2.6 x 10 9 /AH sec -1 cm 3 


(38) 
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where N is the neutral number density and M the average molecular weight. An 
alternate approximation of the ion drag force (Dalgarno, 1964) follows from the 
relation 

v. n /N = 7.3 x 10“ 10 (T/ 1000) 0 • 4 (39) 

Expression (39) yields slightly different results for the ion drag force. A more 
elaborate procedure for the evaluation of ^ is possible, but not necessary as 
the uncertainties of the ion density distribution determine the accuracy. 

In the height region where the collision frequency v in is larger than the 
ion gyro frequency the ion velocity will be close to the neutral velocity and the 
ion drag force will be negligible. This is the case for altitudes below 150 km. 

In the height region where the ion gyro frequency is much larger than the 
collision frequency K n the ions are constrained by the electromagnetic forces. 
There exists a transition region between 150 km to 200 km where the ion motion 
is controlled partially by the neutral motion and partially by the electromagnetic 
forces. In this region the drag forces are difficult to estimate (Lindzen, 1967), 
although the relative low ion densities in this height region tend to decrease the 
importance of the ion drag force in the equations of motion. 

Thus in addition to the uncertainties of the pressure gradients in the lower 
thermosphere the uncertainty of the ion drag force will contribute to the un- 
reliability of the computed wind field in the lower thermosphere. 
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It has been inferred that electric fields in the ionosphere exist due to the 
Sq current distribution. Generally the electric fields will be perpendicular to 
the geomagnetic field and the ion motion perpendicular to the geomagnetic field 
will be given by 

7 ,,„i b/b 2 < 40 > 

The ion motion parallel to the geomagnetic field is determined by the neutral 
motion and given by 

V ionll = (V N • B)5/B2 (41) 

where V N is the neutral air velocity. In the absence of electric fields the ion 
drag force does not depend on the magnitude of B, but only on the direction 
of B. Assuming electric fields do exist and are well-known, then the inclusion 
of the ion drag due to the ion velocity normal to B in our equations of motion 
does not pose any problem. The term [ExB]/B 2 does not involve the neutral 
velocities and one could therefore just add the ion drag due to the electric fields 
to the right hand side of the equations, i.e. to the pressure gradients. 

Electric fields in the ionosphere are difficult to observe directly, so that 
the information about them comes from theoretical deductions of ionospheric 
behaviour. Commonly the electric field of the E-region is extrapolated to 
higher altitudes. From Maeda’s (1971) analysis a lunar component and a solar 
component of the electric potential coexist. Both are of the same order of 
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magnitude but for the steady state wind field only the solar component is re- 
quired. Following Maeda we have calculated the ion velocities normal to B and 
have found them to depend strongly on the higher harmonics of Maeda's analysis. 
For instance if only the P J with 1-2 are included then ion velocities of less 
than 10 m/sec result. If p} with 1-4 are included the ion velocities would 
become as large as 57 m/sec. Inclusion of even high terms of the harmonic 
representation of the electric potentials would probably result in entirely differ- 
ent, and probably larger ion velocities, as the amplitudes of the higher spherical 
harmonics in the representation of the electric potentials do not fall rapidly 
enough to compensate the effects of differentiation. For this reason the ion 
velocities normal to B in the height region of interest to us are not well known. 
Generally (Rishbeth, 1972) it is assumed that they are less than 30 m/sec. In 
the region where the neutral winds are mainly determined by the amplitude of 
the ion drag force, a first approximation of the effect of the ion motion normal 
to B will be simply the addition of the ion velocity to the wind velocity. This is 
seen from the linearized equations of motion. Thus, while the ion motion normal 
to B may be important for particular aspects of the global wind field, like the 
diurnal average zonal motion, i.e. the superrotation of the atmosphere, they are 
not very significant for the general global flow pattern. From these considera- 
tion, especially the uncertainty of the electric fields themselves, we have decided 
not to include the electric fields in our wind computations. With the above 
assumptions the ion drag force f . becomes 
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ff o } n = 2.6 x 10' 9 N . V ((9) //M • cos 2 0/(1 - 0.75 sin 2 0) (42) 

f <*> = 2.6 x 10" 9 N. V^)//M 

i on i 

By dividing the equations of motion (1.4) by u we may define a dimensionless 
ion drag coefficient 


D. on = 2.6 x 10 -9 N. /(/M • o>) . (43) 

The horizontal component of the magnetic field B deviates from the meridional 
direction by the declination angle D. D depends both on longitude and latitude. 

If the dependence of D on longitude is taken into account, then the ion drag forces 
become dependent on universal time. We cannot include in our treatment the 
variations of D without raising the number of independent variables from 3 to 4 
and abandoning the steady state approach. We have therefore to assume that the 
earth’s magnetic dipole is alligned with the earth's axis. This simplification, 
equivalent to setting D = 0 causes the meridional ion drag force to vanish at the 
equator. As a result of the reduced drag in the equatorial zones rather high 
meridional velocities and velocity gradients result from the computation, es- 
pecially at solstice conditions. This is obvious for the linearized equations of 
motion, which become extremely simple at the equator as the Coriolis force 
also vanishes. The non-linear equations avoid this difficulty, but even then the 
meridional velocity gradients are large and the convergence process is slow or 
the iteration scheme may even diverge. Without making assumptions contrary to 
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physical realities we may avoid this difficulty by noting that the meridional ion 
drag depends only on the absolute value of the angle between the magnetic field 
lines and the meridian. We may take the average value of this angle at the equator 
as 11.3°, as a magnetic dipole tilted by 11.3° gives the best fit to the geomagnetic 
field (Mead, 1970). The latitude which corresponds to a magnetic inclination 
I of 11.3° is 5.8°. In line with this reasoning we have substituted at all latitudes 
<E> with | 0 1 < 5.8° in the expression for the meridional ion drag given by (42) 
instead of the factor cos 2 # /(1-0.75 sin 2 #) the value of that factor at a latitude 
of 5.8°. This value is 0.0384. Due to this substitution the meridional ion drag 
force does not vanish at the equator and we have avoided difficulties in the 
convergence process. 

In order to complete the global representation of the ion drag force it is 
necessary to know the ion density as a function Of height, latitude and local time, 
day of the year and solar activity. It will be a challenge for theorist in the future 
to construct a three dimensional model where the ion density is calculated con- 
sistently with the neutral density and the wind system. In our computation we 
have not attempted such a consistent treatment but used a given model of ion 
distribution to determine the drag force. In a early stage of a our computation 
an approximate distribution was used, but in the final calculation we have used 
the Penn State Ionospheric Model which gives numerical values for the ion 
densities as a function of all the above mentioned parameters. In addition, the 
ion density in the Penn State model is dependent on the geographic longitude \. 
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Tliis is not in line with our steady state picture that does not allow functional 
dependence on the longitude alone, but only on the local time ^ - "t. In order to 
adept the Penn State model for our purposes we have averaged the ion densities 
over all longitudes and used these average densities in order to determine the 
drag force. In the actual computation the Fourier coefficients up to second order 
of this longitudinal averaged drag force where used. In this process the longi- 
tudinal averaging reduced the semi-diurnal component of the drag force relative 
to the average drag force to a considerable degree when a comparison was made 
with the drag force that referred to a given fixed longitude. The diurnal component 
of the longitudinal-averaged drag force was also reduced, but somewhat less 
than the semi-diurnal component. The diurnal average of the ion drag coefficient 
is represented in figures 9 and 10 for a solar activity of F 10 7 = 200 at both 
equinox and solstice conditions. Figure 11 shows also the time dependent com- 
ponents of the ion drag for equinox at the equator. The above mentioned figures 
show that the ion drag coefficient, even after our longitudinal averaging, is not 
entirely symmetric at equinoxes. Also the winter and summer solstice coef- 
ficients are not exactly anti-symmetric with respect to the equator. The causes 
of these deviations from symmetry are the various geographical and seasonal 
anomalies of the ion distribution. For these reasons the computed winds will 
also have slight deviations from the symmetry that would be expected from the 
symmetric pressure gradients. 
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The Viscous Forces 

The viscous forces in the equations of motion are given by 

f vis = p/p • (v 2 v + jv(V • V)) (44) 

where p is the dynamic viscosity coefficient and 17 = p/p the kinematic 
viscosity. The term 1/3 V (V . v) is small, it vanishes completely for a constant 
density flow. In neglecting it and also dropping 3 V (r ) /d 9 and B V ( r) /B<£ we 
obtain for the viscous forces in the azimuthal and the meridional direction in 
spherical coordinates (Eskinazi, p. 206): 


f <*> =r> (V 2 V {<i) - 


r (<*>) 


2 cos 9 


r 2 sin 2 9 r 2 sin 2 9 


(45) 


f = -n f V 2 - 


yC 6 *) 


2 cos i9 BV' 


(<*>) 


r 2 s in 2 f? r 2 s in 2 0 ^ 


(46) 


After writing the differential operator explicitly (Eskinazi, p. 207) we may re- 
solve the horizontal viscous force into two parts: vertical viscosity which is 
associated with vertical velocity shears, and horizontal viscosity which is as- 
sociated with horizontal velocity shears. 

The vertical viscosity is approximately equal to the expression 77 B 2 y/B z 2 
which is usually used in an approximate treatment of the equations of motion. 
In our treatment we have included the horizontal viscosity as well (except at 
the poles where it is negligible small). It is given by the expressions 
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if<*) - L_ A 

-n V1S r 2 sin^ 


[s in 9 


be j 


1 ^V (g) V (g) 2 cos 0 BV^_( 47 > 

r 2 sin 2 0 B0 2 r 2 sin 2 0 r 2 sin 2 <9 


if <*> = _L_ i (sin sliiV - 1 . iii . _T!L + ^ 2i!i (48) 

^ r 2 sin(93c? \ b6 ) r 2 sin 2 (9 bcfi 2 r 2 sin 2 0 r 2 sin 2 0 b<fc 

While the vertical viscosity terms are at all latitudes of great importance for 
the resulting flow pattern, the horizontal viscosity terms are only important 
when the meridional or azimuthal velocity gradients are large. Our results 
show that this is indeed the case in a band of latitudes near the equator. The 
order of magnitude of the horizontal viscosity terms becomes in this narrow 
latitude band almost as large as the main terms of the equation of motion. 
Furthermore the inclusion of the horizontal viscosity in the calculation facili- 
ties the convergence of the iteration process, because the horizontal velocity 
gradients are decreased by it and therefore the non-linear terms of the con- 
vective derivative of the velocity become smaller, thereby decreasing the in- 
fluence of the nonlinearities and speeding up the convergence process. 

The kinematic viscosity increases exponentially with altitude because the 
density decreases exponentially. Generally an approximation for the viscosity 
is used (Rishbeth, 1972). In our calculations a somewhat more accurate method 
was applied: The dynamic viscosity coefficient as a function of temperature for 
the various atmospheric constituents was taken from the results of Yun et al. 
(1962) and then the Jacchia model was used to calculate the appropriate tem- 
perature and composition at the latitude, altitude and local time in question. 
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For each altitude and latitude the value of the kinematic viscosity was Fourier- 
analysed and the resulting Fourier coefficients up to second order used in the 
computations. The exact expression for the dynamic viscosity coefficient we 
used was 

< 49 > 

where the summation was for molecular nitrogen, molecular oxygen and atomic 
oxygen, n. are the respective number densities, N the total number density 
and the / j - 0i and a. are 

fx 0 . = 4.017. 10" 4 , 4.771. 10" 4 , 4.771. 10" 4 

a. = 0.62, 0.59, 0.59 i = 1 to 3 

The dependence of the kinematic viscosity in the atmosphere on height is shown 
in figure 12 both for the diurnal average value of the viscosity and for the time- 
dependent components. It is seen that the ratio of the diurnal component to the 
diurnal average increases with height from a value zero at 120 km to about 20% 
at 500 kilometers. This shows the importance of including a diurnal variation 
of the viscosity in the isothermal region. The time of the maximum of the 
viscosity also changes considerably from the lower thermosphere, where it is 
in phase with the temperature, to the isothermal region where it has a phase 
difference of 12 hours relative to the maximum of density. In the upper thermo- 
sphere the diurnal variation of density becomes more important than the diurnal 
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variation of temperature for the evaluation of the kinematic viscosity which is 
approximately proportional to T 0 ' 6 / p . 

Boundary Conditions 

(a) Lower Boundary 

Little is known about the steady state diurnal variation of the atmosphere in 
the 120 km height region. For this reason no observationally founded assumption 
regarding the lower boundary conditions for the global thermospheric wind field 
can be made. The pressure gradients are probably very small in this region, but 
so are the viscous and drag forces. For these reasons a steady state wind velocity 
of the order of 100 m/sec at 120 km cannot be discounted. In our computation we 
have assumed no winds at 120 km in line with the Jacchia model used by us. Fortu- 
nately the effect of a possible wind field at 120 km on the global wind pattern above 
180 km is negligible, as the coupling between adjacent height layers is only through 
the viscosity term of the equation of motion. The kinematic viscosity coefficient 
at 120 km is almost 5 orders of magnitude less than at 500 km (Figure 12), thereby 
decreasing the coupling to a very considerable degree. We performed a test cal- 
culation with a non-vanishing wind field at 120 km. No, or only very little change 
in the resulting wind fields above 160 km were noted as will be shown in Part 2. 
This result is also in accordance with the result of Lindzen (1967). 

(b) Upper Boundary 

The upper boundary conditions are derived generally from considerations 
involving the viscous forces. It is easy to see (Rishbeth, 1972) that 3 2 V/Br 2 -0 
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at exospheric heights is required in order to balance the equations of motion. 

The generally accepted deduction that this implies also ^V/^r - 0 is not as well 
established by strict theoretical considerations. Nevertheless we have also 
used in our calculations the condition 3y/3 r - 0. Ina discussion whether 
boundary conditions with Bv/Br ^ 0 are possible it would be necessary to 
analyse the transition region between the collision controlled height region and 
the collision free region above the exosbase. 

Chapman and Cowling (1952) have discussed the behaviour of the viscosity 
coefficient at low gas densities and have found it to be decreasing under certain 
conditions, but no strict treatment of the transition region at the exosbase re- 
garding the horizontal wind shears is known to the authors, so that a possibility 
of BV/B r / 0 cannot be entirely discounted. In this respect it may be remembered 
that King-Hele (1971) has observed the average azimuthal velocity of the 
atmosphere to be decreasing above 300 km to at least 500 km. Obviously these 
observations cannot be reconciled with a boundary condition By/Br = 0, so that 
observational evidence does also not decisively confirm the assumption BV/3r 
= 0 that is generally made. 

In part 2 the resulting wind field will be discussed. 
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FIGURE CAPTIONS 

Figure 1. Structure of Jacobian Matrix J = BF/BV. according to iteration 

method 1 (basic Newton-Rapphson iteration scheme for one altitude layer 
superimposed Gauss-Seidel iteration over altitude range). For iteration 
method 2 structure would be similar, but there would be 37 blocks of 190 
rows. Each block would correspond to one latitude. Matrix J is tri-diagonal, 
all elements not in blocks indicated are zero. For these elements no com- 
puter storage is required. Sub-matrices typeD k : All elements of these 
matrices correspond to the same altitude. They arise from the various 
terms of the equations of motion. Some terms are also due to viscosity. 
Sub-matrices types L k and U k : These sub-matrices couple adjacent 
altitude layers by the viscosity term r, (B fy/Bz 2 ). They are diagonal in the 
sense that all their sub-blocks of order 10 x 10 are on the diagonal. If 
viscosity would be neglected the matrices and U k would vanish. Sub- 
matrix D ig : This sub-matrix has a special structure due to upper boundary 
condition BV/Bz = 0 at 500 km. 

Figure 2. Fine structure of sub-matrices of type D k shown in Figure 1. Matrices 
D k have sub-blocks of order 10 x 10 denoted by DD k . , DL k . and DU ki . Most 
of the elements of DD kl and DD k37 are zero due to the polar boundary condi- 
tions. The elements of the matrix DD are due to the main terms of the 

k 1 

equation of motion. Some of them are also due to the non-linear terms of 
the convective derivative and the viscosity coupling. Matrices DL and 
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DU k couple various latitudes at a fixed altitude. These sub-matrices would 
vanish when the non-linear terms of the convective derivative and the hori- 
zontal viscosity coupling are neglected. 

Figure 3. Structure of sub-matrices L k and U k of Figure 1. The matrices 

L and U are diagonal in the sense that their sub-blocks, which themselves 

k k 

are matrices of order 10 x 10, are all on the diagonal. The reason for this 
diagonal structure is that viscosity coupling between adjacent altitude layers 
is only considered between mesh points of the same latitude. In method 2 
the equivalent matrices would couple adjacent latitudes at the same altitude 
and would therefore not have this simple structure. Matrices LD ki and 
LD k37 have special form due to boundary conditions at the poles. Matrices 
LD ki couple adjacent height layers due to vertical viscosity. Coupling exists 
only when latitudes are equal. The sub-structure of LD k . shows that BM 
couples the meridional flow and BZ the zonal flow. As the viscous drag 
does not couple meridional to zonal flow when the altitudes differ, the off- 
diagonal sub-blocks of LD vanish. 

k i 

Figure 4. Structure of matrices DL ki , DD k . and DU ki of Figure 2. Computor 
storage for all the 300 elements detailed here is required. Matrices DDM, 
DDZ and DD1 have numbers as their elements. DDM, DDZ and DD1 are 
5x5 matrices, DL is a 10 x 10 matrix. Matrices DD1 couple meridional to 
zonal velocities. Their elements are due to the Coriolis force, the non- 
linear elements of the convective derivative and the horizontal viscosity 
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terms. The off-diagonal elements of matrices DDM and DDZ couple the 
various Fourier modes. If a further linearization of the equations of motion 
is made (additional to the linearization of the convective derivative) and the 
products of higher Fourier modes are neglected then all off-diagonal ele- 
ments of columns 2-5 and rows 2-5 would vanish and the various Fourier 
modes decouple, thus greatly simplifying the computations. 

Figure 5. The height and latitude dependence of the driving force as deduced from 
Jacchia’s model at equinox conditions with a solar activity F ? = 200. The 
diurnal average meridional force and the diurnal amplitudes of the azimuthal 
and meridional forces are shown. The diurnal average of the azimuthal force 
is zero, or nearly zero. The forces are shown for altitudes of 200, 300, 400 
and 500 km. The Jacchia model was modified in the polar region in order 
to fulfill the boundary conditions. 

Figure 6. The same forces as shown in figure 5 for summer solstice conditions. 

Figure 7. The global pattern of the driving force at a height of 300 km for 
equinox conditions. 

Figure 8. The global pattern of the driving forces for summer solstice conditions. 

Figure 9. The latitude dependence of the dimensionless diurnally averaged ion 
drag coefficient for equinox conditions and a solar activity of F 10 7 = 200. 

The coefficients are represented for the altitudes 140, 220, 300, 380 and 
460 km. Ion densities are deduced from the Penn State ionospheric model. 
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Figure 10. The latitude dependence of the diurnally averaged ion drag coefficients 
for summer solstice conditions and a solar activity of F ? = 200. 

Figure 11 . The height dependence of the average and time-dependent Fourier 
coefficients of the ion drag at equinox conditions and a solar activity of 

Fio.7 =20 °- 

Figure 12. The height dependence of the diurnal average and time-dependent 
Fourier coefficients of the kinematic viscosity for equinox conditions at 
the equator as deduced from Jacchia's model at a solar activity of 
F = 200. 
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